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Using the self-consistent tight-binding Bogoliubov-de Gennes formalism we have studied the effect 
of nearest neighbor spin-singlet bond (SB) correlations on Josephson coupling and proximity effect in 
graphene SNS Josephson junctions with conventional s-wave superconducting contacts. For strong 
enough coupling the SB correlations give rise to a superconducting state with either an extended 



2-, or d^y-wave symmetry, or different combinations of the d-waves with the d x i_ i + id x 



state favored in the bulk. Despite the s-wave superconducting state in the contacts, the SB pairing 
state inside the junction has d-wave symmetry and clean, sharp interface junctions resemble a 'bulk- 
meets-bulk' situation with very little interaction between the two different superconducting states. 
In fact, due to a finite-size suppression of the superconducting state, a stronger SB coupling constant 
than in the bulk is needed in order to achieve SB pairing in a junction. For both short clean zigzag 
and armchair junctions the d-wave state that has a zero Josephson coupling to the s-wave state 
is chosen and therefore the Josephson current decreases when a SB pairing state develops in these 
junctions. In more realistic junctions, with smoother doping profiles and atomic scale disorder at 
the interfaces, it is possible to achieve some coupling between the contact s-wave state and the SB 
d-wave states. In addition, by breaking the appropriate lattice symmetry at the interface in order 
to induce the other d-wave state, a nonzero Josephson coupling can be achieved which leads to a 
substantial increase in the Josephson current. We also report on the LDOS of the junctions and 
on a lack of zero energy states at interfaces despite the unconventional order parameters, which we 
attribute to the near degeneracy of the two d-wave solutions and their mixing at a general interface. 

PACS numbers: 74.45.+C, 74.50.+r, 74.20.Mn, 74.20. Rp, 74.70.Wz, 73.20.At 



I. INTRODUCTION 



Since the experimental realization of a single layer of 
graphite, called graphene, a few years agoi the unusual 
two-dimensional Dirac-type spectrum of graphene has 
spurred a lot of attention and has been found to give 
rise to effects not usually expected in condensed mat- 
ter systems^ However, little attention has been given 
to the fact that p7r-bonded planar organic molecules, of 
which graphene is an infinite extension, show a prefer- 
ence for nearest neighbor spin-singlet bonds, or singlet 
bond (SB) in short, over polar configurations, as origi- 
nally captured by Pauling's^ idea of resonating valence 
bonds. In a microscopic Hamiltonian this spin-singlet en- 
hancement takes the form of an intrinsic JSj • Sj term be- 
tween nearest neighbors. Baskaran- combined this term 
with band theory into a phenomenological Hamiltonian 
a few years ago and showed that strong enough coupling 
J will give rise to mean-field superconductivity. The cur- 
rent authors^ extended this work by pointing out that the 
most favorable symmetry for this coupling, and with very 
high mean-field T c , is not, as previously assumed, s-wave 
but instead a two-dimensional <i-wave which in the bulk 
will break time-reversal symmetry. This <i-wave state 
has further been studied within the Dirac Bogoliubov- 
de Gennes (BdG) formalism to calculate properties of 
junctions between d-wave superconducting and normal 



graphene^ Recently many-body approaches have con- 
firmed the possibility of realizing this d-wave state in 
doped graphene. A functional renormalization study 
has found that, for the physical parameters expected in 
graphene, the honeycomb lattice appears to flow toward 
a d + id superconducting stated and variational quantum 
Monte Carlo calculations on the Hubbard model has con- 
firmed the stability of this state toward quantum phase 
fluctuations caused by on-site repulsion^. In addition, by 
using hybrid exchange density functional theory, defect 
arrays in graphene have been shown to cause a ferromag- 
netic ground state even at room temperature for unex- 
pectedly large defect separations. 1 - This state is a con- 
sequence of long-range spin polarization with alternating 
spin direction on adjacent atoms and was interpreted as 
driven by the intrinsic spin-singlet 7r-band instability. 

The mean-field results show that there is a quantum 
critical point at J c = 1.9t at half-filling, and with an 
estimated J ~ t in graphene, no effects of SB corre- 
lations are expected in undoped graphene. But, with 
finite doping any J will give a superconducting state, 
and for enough doping it is expected that T c will be 
measurable. However, the doping levels necessary in 
order to give T c ~ 10 K are beyond the current tech- 
nology of doping graphene by gating. It might be pos- 
sible to achieve a higher degree of doping by chemi- 
cal absorption on graphene, but care has to be taken 
not to alter any other properties of the graphene but 
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the Fermi level. Atomic sulfur absorption on graphene 
seems to fulfill these requirements^ and this offers an 
explanation to the high-temperature superconductivity 
found in graphite-sulfur compositesj 11 ' 12 ' 13 However, as 
a composite this material is not very easily modeled 
and it is therefore logical to seek simpler systems where 
the SB correlations might naturally be enhanced. One 
obvious candidate would be superconducting graphene. 
While a few other proposals exist for superconductiv- 
ity in graphcn o 14 ' 15 ! 16 it has only so far been realized 
inside Josephson junctions. Graphene SNS Josephson 
junctions were recently manufacture d 17 ' 18 ' 19 by deposit- 
ing conventional superconducting metal contacts on top 
of a graphene layer. A finite supercurrent, seemingly in 
agreement with theory^ was measured for these junc- 
tions even at the Dirac point. Our quest here is to in- 
vestigate what possible effects intrinsic SB correlations 
would have in a graphene SNS Josephson junction. More 
specifically, we are interested in the following two ques- 
tions: (1) Is it possible that the superconducting state 
induced from the contacts will cooperate favorably with 
the intrinsic SB correlations to make them stronger and 
therefore easier to detect, and (2) if the SB correlations 
are strong enough to cause superconductivity; what sym- 
metries will be chosen and what are the measurable ef- 
fects, such as changes in the Josephson current and in- 
terfacial zero energy states (ZESs)? 

We will use a self-consistent tight-binding Bogoliubov- 
de Gennes (TB BdG) formulation to address these ques- 
tions. The same framework was recently used by the 
authors^ to study graphene SNS junctions with the SB 
coupling switched off. This study confirmed and comple- 
mented earlier analytical results based on the non-self- 
consistent Dirac BdG formalism^ The benefit of using 
a self-consistent method is that the order parameters are 
allowed to vary spatially. This allows for an explicit cal- 
culation of the proximity-effect depletion and leakage of 
the pairing amplitude around the interface and results in 
a Josephson current properly calculated from the prox- 
imity effect. In addition, when the main interest is the 
possible coupling between different order parameters at 
interfaces as in the current work, a self-consistent method 
is a must. More specifically, we assume the following 
model. The graphene sheet is assumed to be impurity- 
free and ballistic. The influence of the superconducting 
contacts deposited on top is simulated by two effects: 
an on-site attractive Hubbard pairing potential U and 
heavy doping. This will give rise to a conventional, k- 
independent, s-wave superconducting state in the S re- 
gions of the graphene and is the simplest way to simulate 
the effect of the superconducting contacts in the present 
model. Further motivation for this choice can be found 
in Ref. 2l|. The SB correlations are modeled by assuming 
a finite coupling J in the graphene and a finite doping, 
though not nearly as high as in the S regions. Due to the 
significant increase in computational complexity we will 
not solve self-consistently for the chemical potential in 
the system as a function of doping but will assume fixed 



effective chemical potentials fx in the S and N regions. 
We study both clean and smooth SN interfaces as well 
as rough interfaces with atomic scale disorder. Exper- 
iments have indicated a high transparency of graphene 
SN interface s 17 ' 19 so these two types of interfaces should 
in between them capture a realistic situation. Figure 
[T] shows the schematic of (a) the experimental and (b) 
model setup. Note that, even when we are using a large 
enough J to achieve a SB superconducting state in N, we 
will, for simplicity, still call this region N, as in normal 
metal. 
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FIG. 1: Schematic of (a) an experimental SNS graphene 
Josephson junction and (b) the model setup with input pa- 
rameters: on-site pairing potential U , the SB coupling con- 
stant J, effective chemical potential jx, length of normal region 
L, and region where the phase of the order parameter will be 
kept fixed X. 

The paper is organized as follows. In the Sec. II the 
TB-BdG formalism will briefly be explained and the 
choices for the physical parameters will be discussed. 
Next we will report our results. First, it is important 
to understand the effect of the two coupling parameters, 
U and J, in the translational invariant bulk, and how 
their superconducting states can influence each other by 
induced pair amplitudes as well as Josephson coupling 
between them at an interface. We then show the re- 
sults for clean zigzag and armchair interfaces. Following 
these results, we expand the study to extended inter- 
faces, including smoother doping profiles and interface 
roughness, and show that additional effects can occur. 
We then discuss local density of states (LDOS) results 
for the junctions and the presence of ZESs. In the last 
section we summarize our results and discuss possible ex- 
tensions and experimental opportunities. 
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II. METHOD 
A. Tight-binding BdG 

The method used here follows closely that of our ear- 
lier workSi where graphene SNS Josephson junctions 
were studied but the spin-singlet bond correlations in 
graphene switched off. In fact, the only difference to the 
formalism is the addition of the SB pairing term but for 
self-containment we choose to include a brief summary 
of the full formalism. 

Following the motivation in Sec. I, we model the SNS 
Josephson junction using the following Hamiltonian: 



t E UU* + alfia) + E HWLfa + gin*) 

<i,j>,o- i.cr 

E ni: ^/Wi>. • </r//i •</!>.• 

i 

E 2J(i)4'Hl (!) 



<i,j> 



Here /j (g ia ) is the creation operator on the A (B)-site 
in cell i of the honeycomb lattice, see Fig. [2] The first 




FIG. 2: (Color online) The graphene honeycomb lattice with 
the two different atomic sites / and g, the three nearest neigh- 
bor vectors {ai,a2,a3} (marking bonds 1, 2, and 3, respec- 
tively), and the zigzag and armchair interfaces marked. 

two terms give the tight-binding band structure for the 
two 7r-bands with hopping parameter t — 2.5 eV and 
a site-dependent effective chemical potential /2(i) which 
will determine the local occupancy, jl = corresponds 
to the Dirac point. The third term models the influence 
of the superconducting metal contacts on the underlying 
graphene by an on-site attractive Hubbard £/-tcrm. This 
term will only be nonzero in the S regions. Finally, the 
last term describes the SB correlations. Here 
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V2 



(/iWi-/Wt) 



(2) 



is the spin-singlet creation operator. We have for sim- 
plicity defined the strength of the J-term as given by its 
value on the /-atom in each spin-singlet. The factor of 



2 comes from the two atoms per unit cell. As will be 
motivated below, J is only set to be nonzero in the N 
regions. The spin-singlet coupling can also be written as 



Jh\.hii = J ( Si-Sj - ~rnn } J , 



(3) 

which shows the similarity between our phcnomcnolog- 
ical model and a t-J model, though in the latter case 
double occupancy is explicitly forbidden, whereas here 
we do not impose any such restrictions. 

To proceed, we use the Hartree-Fock-Bogoliubov 
mean-field approximation and arrive at the following 
quadratic Hamiltonian: 

H M F = ~ * E (/iUi+a,ff + aLfi-a,*) 
<i,a>,(T 

+E^ i )^/-+9 i W) 

i,<7 

+ EA C /(i)(/ i t T / i t i + .g i t T5 / i ) + H.c. 

+ E A .'a(/iW+al - /iW+a T ) + H.C., (4) 
i.a 

where we have introduced the three nearest neighbor vec- 
tors an, with i = 1, 2, 3, to index the nearest neighbor di- 
rections, see Fig. [2j The spatially dependent mean-field 
superconducting order parameters are defined as 



A v (i) = -U(i) 



(MM) + (fiifit) 



Aja(i) = -J(i)(/i;£?i + aT - fiWi 



+al/- 



(5) 



The order parameter for the SB coupling is depen- 
dent on a.i and we will treat all three nearest neigh- 
bor bond directions independently in order to include 
all possible symmetries of this order parameter. We 
will use A ,/ = ^|A, 7ai | 2 + |A Ja2 | 2 + |A Ja3 | 2 as a mea- 
sure of the strength of the SB order parameter. In or- 
der to study proximity effects we also need the pair- 
ing amplitudes for each order parameter which are sim- 
ply defined as Fu(i) = -Au(i)/U(i) and F J(a) (i) = 

-A JW (i)/(V2J(i)). 

The mean-field Hamiltonian can be rewritten in the 
standard BdG formalism (see Refs. [HHHl] for recent 
applications) which in a tight-binding formulation con- 
sists of an eigenvalue problem which needs to be solved 
for its positive eigenvalues v. 



E 



( ffo(ij) A(i,j) 
^At(i,j) -Fo(iJ) 



= E p 



(6) 



For the two-atom graphene unit cell the matrices H and 
A will have dimensions 2x2 and can be written as 



ffo(i,j) 
A(i,j) 



M(i)<*U 



E a Aja<5i- 



-* Ea ^i+aj 
aj HWii 

A u(i)Sij Ea A ./a^i+a,j 
a.j 



A[/(i)<5ij 



(7) 

(8) 
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Thus N unit cells leads to a AN x AN eigenvalue problem 
of which 27V eigenvalues will be positive since H = 
and A = A T . Writing the eigenvectors explicitly as 
four-dimensional (per unit cell i) using the new nota- 
tion (u, v) —>■ (u, y, v, z) we can write the self-consistency 
equations as 

A u(i) = ^f- + Vi*n tanh ^ (9) 

A Ja (i) = J(i) X)(wr+a»r + tanh ^— . (10) 

Any structure with spatially varying //, U, and J can 
now be solved by starting with an initial guess for Ajj and 
A,/ a . After finding the 2N positive eigenvalues of Eq. (J6j) 
we can compute new values for Ajj and Aj a using Eqs. 
© and (p~0|) and repeat the process until the differences 
in Ajj and Aj a between two subsequent iterations are 
less than the desired accuracy and thus self-consistency 
is reached. 

From a computational point of view it is beneficial to 
reduce the size of the eigenvalue problem as much as 
possible. This can be achieved by assuming translational 
symmetry in the direction perpendicular to the junction 
and then applying Bloch's theorem. We have imple- 
mented both clean, smooth interfaces, with full transla- 
tional symmetry along the junction interfaces, and simu- 
lated the effect of a rough interface by using an extended 
unit cell which includes a total of four graphene unit cells. 

One of the main characteristics of a SNS junction is 
the Josephson current it can carry when a finite phase 
gradient is applied across it. Here, since the contacts 
are modeled with the attractive [/-term, a phase differ- 
ence needs to be established in Ajj in order to produce 
any supercurrent. We achieve this by fixing the phase of 
Ajj in the outermost regions, labeled X in Fig. [I] of the 
contacts and then solving self-consistently for both phase 
and amplitude in the rest of the structure. The current 
can then be calculated relatively straight-forward using 
the continuity equation and the Heisenberg equation for 
the electron density, see Ref. [2l| for further details. 

The main output information from a self-consistent so- 
lution of the TB BdG formalism is (1) the pairing am- 
plitudes Fjj and Fj a which will describe the proximity 
effect in the junction and (2) the Josephson current the 
junction carries when a finite phase gradient is applied 
across the junction. 

B. Simulation details 

With the directional dependence of the SB coupling, 
different interfacial directions can exhibit different behav- 
iors. We study the two main interfaces for the graphene 
honeycomb lattice which are the zigzag and the armchair 
terminations, see Fig.O Interfaces in a random direction 
are called chiral interfaces and while we do not explicitly 



study these we will be able to infer general results from 
the two limiting cases of zigzag and armchair interfaces. 

The physical input parameters are the on-site pairing 
potential U in S, the SB pairing potential J, the effective 
potential jl in S and N, the length L of N, and temper- 
ature. For the superconducting contacts we choose the 
following setup: U(S) = 3.4 eV = 1.36*, fl(S) = 1.5 eV = 
0.6t. This leads to Ajj — 0.1 eV and a superconducting 
coherence length £ ss 50 A which corresponds to 25 unit 
cells in the zigzag direction and 40 unit cells in the arm- 
chair direction. These values satisfy Xp(S) <C £, allow 
us to numerically investigate both the L < £ and L > £ 
cases, and coincide with the choice we used previously 21 
when investigating SNS junctions without the J-term. 
They are however quite large values for a realistic sit- 
uation but smaller superconducting gaps lead to slower 
convergence rates and also the need for larger systems 
making calculations less feasible. We have checked our 
key results for smaller U and found no significant differ- 
ence. 

We choose to set J = in the S regions. Physically 
this can be motivated by the fact that the graphene is 
here in very close proximity to a large metallic electron 
reservoir which provide ample screening of the on-site re- 
pulsion which in turn leads to significantly reduced effec- 
tive coupling J. This setup is not akin to the situation 
in the intercalated graphites where we previously have 
argued^ for diminishing SB correlations. Also, as we will 
discuss below, from a practical point of view, any rea- 
sonably value of J together with a large attractive U will 
not change the physical state in any other way than to 
induce an even larger Ajj. We initially studied junctions 
with the physically realistic valued J = t = 2.5 eV and 
with doping levels in N ranging from zero to no Fermi 
level mismatch (FLM) at the interfaces, i.e. /i(N) = p,(S). 
This value of J gives T c ss 10 K for fi = 0.7 eV in the 
bulk. However, as will be clear below, J = t only induces 
a superconducting state in either very long, approaching 
bulk conditions, or in very heavily doped (Jl > 0.7 eV) 
junctions. In order to keep the system sizes down and the 
doping level such that a significant FLM is still present 
at the interface, but still study the effect of order param- 
eter symmetry selections and the physical consequences 
thereof, we have also preformed simulations with J « 1.5t 
and /z(N) = 0.7 eV. Different doping levels in S and N 
are especially important in order to study effects on the 
Josephson current as no FLM in fact leaves the whole 
junction fully superconducting at the currently accessi- 
ble junction lengths^ 

For most junctions we have for simplicity assumed that 
the doping profile changes abruptly from p,(S) to /i(N) at 
the interface. Physically this means insignificant doping 
leakage from the S regions to the N region. To exam- 
ine the effect of this approximation we have also stud- 
ied junctions with a linear doping profile drop over a 
few unit cells. In terms of L, we have studied zigzag 
junctions with L = 10, 30, 60 unit cells (1 unit cell = 
2.13 A) and armchair junctions with the corresponding 
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L = 18,52,104 unit cells (1 unit cell = 1.23 A). The 
temperature was chosen to be T = 10 K throughout the 
work, which in comparison to T c in the S regions is effec- 
tively zero temperature. 

The accuracy of the solution is determined by the 
choice of termination criterion for the self-consistency 
step, the number of Appoints used in the Fourier trans- 
form preformed in the direction parallel to the junction, 
and the size of S and X to ensure bulklike superconduct- 
ing conditions in the S regions, all of which have been 
tested thoroughly. 

III. RESULTS 

Before proceeding with the numerical results for SNS 
Josephson junctions it is important to understand what 
different symmetries the order parameter for SB pairing 
can take in the bulk and also how these can couple to the 
s-wave symmetric superconducting state induced from 
the contacts. 



A. Order parameters in the bulk 

In the clean limit, the bulk has translational symmetry 
in both spatial directions. This reduces Eq. ^ to a 1 x 1 
matrix eigenvalue problem, which needs to be solved for 
each reciprocal vector. We will here briefly review the 
different superconducting states possible when U and/or 
J 7^ 0, results that have all been published elsewhere, as 
well as how nonzero pairing amplitudes can be induced 
through order parameter interaction. 

1. [7/0, J = 

This case corresponds to the attractive Hubbard model 
on the graphene honeycomb lattice which mean-field so- 
lution is straightforward to obtain (see e.g. Ref. [Ill)- At 
half-filling, i.e. fl = 0, the zero density of states leads 
to a quantum critical point such that U > U c = 2M is 
necessary for a superconducting state. However, at fi- 
nite doping any U will give a finite superconducting T c . 
Since the pairing is on-site, no fc-dependence is present 
and the order parameter has a conventional s-wave sym- 
metry throughout the Brillouin zone. 

2. U = 0,J^0 

/ / represents the case of SB correlations which 
are present in all p7r-bonded planar organic molecules. 
By treating the three nearest neighbor bonds indepen- 
dently there are three degrees of freedom for the super- 
conducting order parameter Aj a . The mean-field solu- 
tion of this model has been studied by the authors in 
a previous paper— and it was found that the solution 



can either have an extended s-wave symmetry or be- 
long to a two-dimensional d-wave symmetry class. The 
superconducting order parameter in graphene will be- 
long to the Dq classi with the s-wave belonging to the 
A\ irreducible representation and the d-waves to the 
Ei representation. The extended s-wave is found when 
A Ja = (A Jai , Aj a2 , Aja 3 ) oc (1,1,1), i.e. when all three 
bonds have the same weight. It is an extended s-wave in 
the sense that the order parameter has the same symme- 
try as the band structure, | ^ a e* k ' a |, in the band struc- 
ture Brillouin zone, i.e. in the reciprocal space where 
the kinetic energy term is diagonal. This s-wave has 
only nodes at the corners of the Brillouin zone. The 
d-wave solution is two-dimensional and is spanned by 
{(2,-1,-1), (0, 1, -1)}. Here the (2, -1, -1) choice has 
a d x 2_ y 2 symmetry in the band structure Brillouin zone 
whereas the (0,1,-1) has a d xy symmetry, see Fig. [3] 
For simplicity we will use the notation d\ = d x 2_ v 2 and 
d-2 = dxy At T c they are degenerate, but below T c the 
complex combination d\ + id2 has the lowest free energy, 
which means that the superconducting state breaks time- 
reversal symmetry. At half-filling the s-wave and d-wave 
solutions are degenerate and, as in the case of the attrac- 
tive Hubbard model, there exists a quantum critical point 
at J c = 1.9t such that J > J c is necessary for a mean- 
field superconducting state. This is significantly higher 
than the estimated J ~ t. However, for finite doping 
any J will give a superconducting state and the d-wave 
solutions will be favored for any reasonable values of J 
and jl. A typical value of the mean-field T c for J = t, 
fx = 0.7 eV is 10 K. 



3. U / 0, J / 

When both U and J are nonzero the order parameters 
Ajj and Aj a have the possibility of mixing. Since we are 
considering U ^ only in the S regions where also the 
doping level is very high, we can simplify the Hamilto- 
nian in this case to only include one of the 7r-bands. It is 
then straightforward to show that among the four possi- 
ble order parameter states the two d-wave states from the 
J-term are left unchanged by a finite U and also Aj/ = 
for these states. The remaining two states are mixtures 
of the two s-wave states, i.e. Ay ^ and Aj a oc (1, 1, 1). 
In terms of T c the mixed s-waves states are heavily ben- 
efiting from a finite U . In fact, for J = t and strong 
enough U to create a stable induced s-wave from the 
contacts, one of the mixed s-wave states has a higher T c 
than the d-wave states. This leads to a situation in the S 
regions where, from a simulations point of view, having 
J zero or nonzero is irrelevant, especially since, as seen in 
Sec. IIIA4, a finite attractive U will always induce a finite 
s-wave pairing amplitude Fj a so even proximity-type ef- 
fects are not ignored by setting J = 0. For simplicity we 
therefore, unless otherwise stated, never assign both U 
and J nonzero at the same site. 



6 



4- Induced pairing amplitudes 

In the above sections we discussed the different or- 
der parameters possible when either U or J are nonzero. 
When studying the proximity effect we are, however, also 
interested in induced pairing amplitudes from one order 
parameter to the other. For U ^ it is straightforward to 
show that a finite Fj a with s-wave symmetry will always 
be produced. Thus inside the S regions where the super- 
conducting contacts produce an on-site s-wave state, a 
SB s-wave amplitude will also be present. Conversely, a 
finite SB s-wave state will induce a finite on-site pairing 
amplitude. However, a finite SB d-wave state has a zero 
overlap with the on-site s-wave state and, therefore, as 
soon as a (i-wave state is chosen when J ^ 0, there is 
no effect on the on-site pairing amplitude in the bulk. 
This can easily be understood from the simple fact that 
the sum over the whole Brillouin zone of a function with 
sixfold symmetry multiplied with a function of fourfold 
symmetry is zero. 



can instead be written as J\ sin(A0) + J2 sin(A0 + n/2) 
where the terms Ji depend on the symmetry of the d\- 
and d2-wave components, respectively, at the interface. 

For our purpose it is enough to understand the possi- 
ble couplings in a SS' junction on a schematic level only. 
Figure [3] shows the three different cases; an on-site s- 
wave on the left-hand side and the three different SB 
order parameters on the right-hand side. As for any elec- 
tronic property it is most important to consider the situ- 
ation around the Fermi surface, which in doped graphene 
consists of two nonequivalent circles at two neighboring 
corners of the Brillouin zone. Figure [3] explicitly shows 
that while the d±- and the e^-wave states have a four-fold 
symmetry in the whole Brillouin zone they have effective 
p y - and p^-wave symmetries on the Fermi surface, re- 
spectively. Note that this is the same <i-wave state as 
discussed in Ref. 0] but there the low energy proper- 
ties are instead given when using atomic operators. We 
have instead used the band structure basis where the ki- 
netic energy is diagonal. As can be seen in the figure, 



B. Order parameter coupling at interfaces 

So far we have only discussed what happens in the 
translational invariant bulk but in order to interpret the 
SNS junction results it is also important to understand 
how two bulk order parameters would couple to each 
other if lined up at an interface. The most relevant sit- 
uation for the SNS structures studied here is obviously 
the possible coupling between an on-site, fc-independent, 
s-wave state on one side and a SB s- or e?-wave state on 
the other side of the interface. 

This situation corresponds to the Josephson effect in a 
SS' junction. While we do not have an explicit insulator 
or barrier in between the two superconductors, the FLM 
at the interface will act as an effective barrier, as will be 
clear from the results below. Quite generally, if we write 
the order parameters on each side of the junction as 

/ ftfe^*, forS' 
I n£e" p k , for S, 

where k = 1,2 allows for two component order param- 
eters, then the supercurrent density can be calculated 
using Ginzburg-Landau theory a a 25 > 26 

2 

J= J2 J ckiM0k -0k), (12) 
fe,i=i 

where J c ki oc hf n^. The consequences in our case de- 
pend on the dimensionality of the SB order parameter. 
For the extended s-wave or a one component rf-wave 
state, the magnitude of the coupling depends on the sym- 
metry of that state at the interface and the current-phase 
relation has the usual sin(A<^)) form, where A<f> is the rel- 
ative phase of the two superconductors. For the time- 
reversal symmetry breaking d% + icfe-wave the current 
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FIG. 3: (Color online) Order parameter symmetry in the Bril- 
louin zone and projected to the Fermi surface at finite doping; 
positive order parameter (blue), negative (red). On the S side 
(left) the order parameter is conventional s-wave but on the N 
side (right) for strong enough J we can have either extended 
s-, di-, (fe-wave or any complex combination of the d-waves. 
The directions of transport for zigzag and armchair interfaces 
are indicated with large arrows. Note that the schematic does 
not include a possible additional overall phase difference be- 
tween the S and N sides. 

for transport through a zigzag interface, a partial over- 
lap of the order parameters on the Fermi surface exists 
for the s- and ^-states as long as we exclude intervalley 
scattering but the overlap is zero for the <ii-state. The 
opposite relation is true for the <i-wave states in trans- 
port through an armchair interface. We would therefore 
expect the Josephson supercurrent through a zigzag SNS 
junction to be enhanced if the N region is in a d%- or 
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extended s-wave SB superconducting state or similarly 
for an armchair junction if the N region has d\- or ex- 
tended s-wave symmetry, but that no effect will be seen 
for the other choices of ci-wave symmetry. For the two 
component d\ + idi order parameter there is finite over- 
lap only for the di-part at an armchair interface and only 
for the d2-part at a zigzag interface. However, note that 
in this particular case, due to symmetry, the extra tt/2 
phase shift will effectively be part of the overall phase 
difference. 

It might be worth noting before closing this sec- 
tion, that the above graphical representation based on 
Ginzburg-Landau theory only treats the coupling be- 
tween the two superconducting order parameters to first 
order. It has been shown by treating the tunneling pro- 
cess to all orders in perturbation theory that even if the 
first-order term, proportional to sin(A0), disappears for 
a SS' junction there might still be contributions from 
higher order terms, starting with a sin(2A0) term.— 



C. SNS junctions with clean interfaces 

With the preparation of the above sections, the inter- 
pretation of the numerical results for SNS junctions is 
quite straightforward. We start with examining clean 
interfaces where the doping profile drops sharply from 
/i(S) to /x(N) at the interface and the ?7-term is abruptly 
interchanged for a J-term. 

Figure 0] shows typical results for a doped zigzag junc- 
tion for different lengths L. Here Aj, along with its 
character, in the N region in a SNS junction is shown 
together with Aj for a slab of the same length. The slab 
has been left with dangling bonds at the surfaces. As 
can be seen, despite the presence of the Ajj supercon- 
ducting state in the contacts, Aj is not enhanced in an 
SNS junction over the value in the slab. In addition, Aj 
is always suppressed in finite-size zigzag slabs compared 
to the bulk value. In terms of symmetries, we see that 
even the smallest slabs have a (ii-wave symmetry. While 
we would expect a d± + ic?2-wave symmetry in the bulk, 
the preference for the di-wave for finite-size zigzag slabs 
can be understood when considering that bonds 2 and 3 
are equivalent whereas bond 1 is not for this interface. 
Since d\ cx (2,-1,-1) it becomes energetically favored 
over the G?2-wave solution where bonds 2 and 3 instead 
have a 7r-phase difference. In fact, this surface/interface 
effect last for the longest zigzag slabs (~ 500 A) we can 
model. The only noticeable effect of the superconducting 
contacts in a zigzag junction is the s-wave symmetry of 
Aj for weak SB pairing. However, as seen in the inset of 
Fig.[D^b), the s-wave state is always too weak to influence 
the Josephson current through the junction. This weak 
extended s-wave state appears only when J and ft are too 
weak to cause a notable (d-wave) superconducting state 
in N and in fact only for finite doping levels. In undoped 
or barely doped graphene, Aj will always have d-wave 
character for any J (including J = 0), but of course, the 
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FIG. 4: (Color online) (a) Aj, (b) s-wave, and (c) di-wave 
character of Aj for zigzag slabs (dashed) and SNS junctions 
(solid) for L — 10 (black), 30 (green), and 60 (red) unit cells 
when J = 4 eV and p,(N) = 0.7 eV. The bulk value of Aj for 
the di-wave symmetry is indicated with a dotted line. 



amplitude of the order parameter is (vanishingly) small. 
So despite the possible coupling between the S region 
on-site s-wave and the extended s-wave state of the SB 
pairing, the lowest energy state when SB pairing devel- 
ops in the junction is still a d-wave state, as in the bulk. 
This applies even for very small junctions. 

For the armchair interface Fig.[5]shows the correspond- 
ing picture. There are a few differences compared to the 
zigzag interface. First, while even the longest slab has d\- 
wave symmetry at the surfaces, it goes toward d\ + id%- 
wave symmetry (di + idi — 50% di-wave character and 
50% c?2-wave character) in the interior. This can be un- 
derstood when examining the armchair interface more 
closely. Here, bond 2 in cell i is only equivalent to bond 
3 in cell (i + 1) and only if the paring potential param- 
eters U and J are not changed in between the two cells. 
This is a much weaker symmetry condition than for the 
zigzag interface. Second, Aj is enhanced in finite-sized 
slabs compared to the bulk. This is because the value of 
A,/ ai is very high at the surface. We have checked that 
this is not due to a localized surface state on this bond. 
If one was to suppress this value at the surface, Aj for 
a finite-sized slab will again resemble the value in a SNS 
junction. Since the value of Aj ai is not allowed to be as 
high at the interface with the S region as it is at a dan- 
gling bond surface the size of Aj in the SNS junction can 
again be understood in terms of the finite-size suppres- 
sion in a simple slab. The symmetry choice for the slab 
is independent of any suppression of Aj ai at the surface. 
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FIG. 5: (Color online) (a) Aj, (b) s-wave, and (c) di-wave 
character of Aj for armchair slabs (dashed) and SNS junc- 
tions (solid) for L — 18 (black), 52 (green), and 104 (red) unit 
cells when J = 4 eV and fl(N) = 0.7 eV. The bulk value of 
A j for the di + id2-wave symmetry is indicated with a dotted 
line. 



Third, the s-wave state lasts for somewhat longer junc- 
tions, or, alternatively larger J-values, in an armchair 
SNS junction, though its effect is similarly negligible on 
the Josephson coupling due to its smallness. Fourth, for 
intermediate lengths L the armchair SNS junction shows 
c^-wave symmetry while for even longer junctions (not 
shown in figure) there is a d\ -Md2-state in all of N. Thus 
the only difference between an armchair slab and a SNS 
junction in terms of Aj is the extended s-wave in short 
junctions and the choice of dy- or d\ + zc^-wave symme- 
tries when the slab instead shows di-wave in the whole 
junction or only at the surface, respectively. 

Figures |4] and [5] establish the close connection between 
a simple slab and the N region in a SNS junction for 
the SB order parameter but we are also interested in the 
development of the on-site pairing amplitude and the su- 
percurrent through a SNS junction with varying J. Fig- 
ures inja) and [H^b) show how these properties develop 
with increased J in a zigzag SNS junction. We clearly 
see that for all </, causing either extended s-wave or d\- 
wave symmetries in N, there is no effect on the amplitude 
of the on-site pairing amplitude Fjj] its proximity effect 
depletion in S and leakage into N is unaffected by any SB 
correlations or superconducting state in N. Interestingly, 
the proximity length of the SB pairing as displayed in 
Fig. [B](c) is quite different from that of the on-site pair- 
ing. The SB pairing amplitude rises very abruptly at the 
interface and no effect of J is seen in the S regions. As 
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FIG. 6: (Color online) Zigzag SNS junction with L = 10 
unit cells for four different values of J in N: J = (black), 
3.5 (green), 3.75 (red), 4 eV (blue) and with fi(N) = 0.7 eV. 
Absolute value (a) and phase (b) of the on-site pair amplitude 
Fu, absolute value of the SB pair amplitude Aj (c), and d\- 
wave character of Fj (d). Remaining character is s-wave. 
Vertical lines indicated the interfaces. Inset in (b) shows how 
the current changes with J relative to the current at J = 0. 



seen in Fig. [H the proximity length for the SB pairing is 
quite short at a surface as well, but the difference is even 
more pronounced in the SNS structure. When J is strong 
enough to cause a di-wave state in N, we see a clear sig- 
nature on the phase of Fjj. The total phase drop over 
the junction remains the same, and thus the / vs A(j) re- 
lationship is unchanged, but there is a significant 'bulge' 
developing in N for aig(Fu). The inset in Fig.[5fb) shows 
the development of the Josephson current with increased 
J compared to when J = 0. For all s-wave states there 
is a small increase in the current (less than 1%) caused 
by the coupling between the extended s-wave in N and 
Ay in S. But since Aj is very small in this situation the 
increase in current is also quite small. Once the di-w&ve 
state develops we see a significant decrease in the current. 
This can be explained by the following argument: the 
e?i-wave state has according to Fig. [3] zero overlap with 
a s-wave state in a SS' junction, so there is no coupling 
between the two condensates, at least not to first or- 
der. Excluding possible contributions from higher order 
terms in this coupling, the current carried through the 
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junction can only have two origins; quasiparticle tunnel- 
ing and Andreev tunneling in which an incoming Cooper 
pair on the interface in S is transformed into an out- 
going electron and an incoming hole on the interface in 
N. These are the same physical processes present in a 
regular SNS junction and the latter process should dom- 
inate at the low temperatures we are considering. Both 
of these tunneling processes are sensitive to changes in 
the normal or, in the case of a superconducting state, 
quasiparticle DOS in N. For J = the N region is just 
doped graphene which has a metallic DOS. When the 
di-wave state develops, the DOS is changed to a nodal 
spectrum since the di-wave has nodes on the Fermi sur- 
face. This effectively decreases the available DOS for the 
above mentioned tunneling processes and the current will 
accordingly decrease. The same mechanism will decrease 
the tunneling current in the case of a s-wave state in N 
as well, but in this case the superconducting energy gap 
is very small and the decrease is also counteracted by the 
Josephson coupling between the two different s-wave or- 
der parameters. The final result is a marginal increase in 
the current instead. 

The corresponding picture for an armchair junction is 
very similar, except for the symmetry of the <i-wave state. 
As commented above, for intermediate length junctions 
and strong enough J, a c^-wave state will develop in 
the armchair junction. This state has zero overlap with 
the S region s-wave state and we again see a drop in 
the current as for the zigzag junction. For the longest 
junctions we could study the zigzag junction is still in a 
di-wave state whereas the armchair has a d\ + id,2-~w&ve 
symmetry. This state has a finite overlap of the <ii-wave 
part with the S region s-wave state but counteracting 
this increase in the current is the fact that the d\ + id,2- 
wave state has a full gap on the Fermi surface and, for 
the J-values we are considering, also a considerable size 
of the gap. This effectively diminishes any quasiparticle 
and Andreev contribution to the current. The final result 
is in fact a reasonably large decrease in the total current 
compared to when J — 0. For the armchair interface the 
'bulge' structure on the phase of Fu is present for both 
the e?2-wave and the d\ + ic?2-wave order parameters. 

The choice of symmetry for the SB order parameter in 
N is interesting. For zigzag junctions, clearly the surface 
preference of bond 2 = bond 3 ^ bond 1 will dictate a 
di-wave state for even quite long junctions. On the other 
hand for armchair junctions, a slab will always have a d\- 
wave symmetry at the surface but in a SNS junction the 
symmetry is either d^- or d\ + ie?2~wave in the whole N 
region. It is almost as if the system has the lowest free 
energy when it can, to as large extent as possible, avoid 
overlap between the S region s-wave and the N region 
d-wave order parameters. 

From the results in Figs. HUB] we can draw the follow- 
ing three conclusions of the effect of SB correlations in 
clean graphene SNS Josephson junctions: (1) SB pairing 
in a SNS junction is going to suffer from finite-size sup- 
pression, in which a larger J (or p) than in the bulk is 



necessary in order to achieve any SB pairing. This clearly 
answers our first question as to whether a conventional 
SNS junction can enhance the effect of SB correlations or 
not. (2) Even with a large enough J such that a SB su- 
perconducting state develops in N, the junction is going 
to resemble the situation 'bulk-meets-bulk' in terms of 
the lack of induced effects in between the different order 
parameters. (3) The current through the junction will 
depend on both the possible overlap of the different or- 
der parameters in S and N but also to a significant degree 
on the DOS of quasiparticles in N. 

Since these conclusions, with the exception of the pos- 
sible symmetries for the d-wave state in N, are identical 
for the zigzag and the armchair interfaces we expect that 
they will be true for any clean chiral interface SNS junc- 
tion. The issue of choice for the d-wave state will require 
a detailed knowledge of the interface structure. 



D. SNS junctions with rough interfaces 

After studying the clean interface SNS Josephson junc- 
tions we are still left with the question if it is possible to 
induce a coupling between the on-site s-wave Ajj and 
the SB (i-wave Aj by making the interface less idealis- 
tic. We first investigate the effect of the sharp doping 
level drop at the interface. It is reasonable to assume 
that some doping is going to leak into N from the higher 
doped S region causing either the change /Et(S) <-> /x(N) to 
take place at another site than where [/^O^J^Oor, 
even more realistically, assume a steady, nonabrupt, drop 
from /2(S) to /2(N) over a few unit cells. We have inves- 
tigated both possibilities and they show similar features. 
Figure [7] shows the situation where the chemical poten- 
tial is allowed to drop linearly from the value in S to the 
value in N over 10 unit cells. As seen, a finite Aj gives 
rise to some induced on-site pairing amplitude Fjj . Note 
however that the much larger value of the SB pairing 
amplitude Fj in the region where the doping drops take 
place is at least primarily due to the effectively higher 
chemical potential in that region and not because of cou- 
pling to Fjj. A similar effect can be seen even with a 
sharp drop in chemical potential as long as it does not 
coincide with the change from on-site to SB pairing po- 
tential. This behavior is different from the clean interface 
and indicates that the big FLM at the interface is a key 
component in the suppression of coupling between the 
two order parameters. But, even with the removal of 
this constraint, there is still a finite-size effect present, 
i.e. a larger J than in the bulk is still needed in order 
to induce a finite Aj and hence to see an increased F(j. 
Also, the symmetry choice for Aj does not change and 
the intriguing 'bulge' structure is still present. In terms 
of the Josephson current through the junction, it is obvi- 
ously increased when Fu is increased, in this particular 
case by about 10% compared to when J = 0. 

Figure [7] shows that it is possible to get induced pairing 
amplitude effects between the on-site s-wave state and a 
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FIG. 7: (Color online) Zigzag SNS junction with linear drop 
in doping profile over 10 unit cells from p,(S) = 1.5 eV to 
p.(N) = 0.7 eV for J = (black), 4 eV (red). Absolute value 
(a) and phase (b) of the on-site pair amplitude Fu, absolute 
value of the SB pair amplitude Fj (c), and di-wave character 
of A j (d). Remaining character is s-wave. Vertical lines 
indicate the SN interfaces as well as where the linear doping 
profile drop ends. 



SB d-wave state. An additional question is if it, by inter- 
face engineering, is possible to choose the d-wave state in 
N that has a nonzero overlap with the s-wave state at the 
interface, and thereby allow for coupling between the two 
supercurrents. For the zigzag interface this would mean 
developing a d2-wave symmetry in N. Figure [8] shows that 
this is indeed possible when care is taken in the interface 
region. More specifically, the strong symmetry between 
bond 2 and 3 has to be broken. We achieve this by dou- 
bling the unit cell in the direction along the interfaces and 
setting U and J to different values on the two /-atoms 
in this new extended cell. Here the green (light grey) 
curve represents a junction where the c^-wave symmetry 
is chosen and the current is also increased by 50%. There 
is also a reasonable inducement of on-site pairing ampli- 
tude from the SB pairing. In order to establish that the 
current increase is not just due to this increase in on-site 
pairing but also due to coupling of the two condensates, 
we compare the results with a system where J 7^ is 
allowed to overlap with U 7^ 0. In this system, plotted 
in red (grey), there is a very strong coupling between the 



FIG. 8: (Color online) Rough interfaces for zigzag SNS junc- 
tion for J = (black), J = 4 eV (red, green) with a sharp 
doping profile drop to /i(N) = 0.7 eV at the innermost ver- 
tical lines. In the interface regions (between vertical lines at 
each interface) U 7^ for all cells and J ^ for 50% of the 
cells (red) and U = 0, J 7^ or U 7^ 0, J = in equal mixture 
(green). Absolute value (a) and phase (b) of the on-site pair 
amplitude Fu, absolute value of the SB pair amplitude Fj 

(c) , and d-wave character, d\ (solid), c?2-wave (dashed) of Aj 

(d) . Remaining character is s-wave. 



on-site s-wave and a SB s-wave in the interface region. 
Despite the much larger induced Fu in this latter case 
the current increase is only around 20%. This proves 
that the majority of the current increase in the d,2-~w&ve 
case comes from the overlap of order parameters and the 
flow of supercurrent between them. 

Interestingly we also see that the 'bulge' on arg(Fu) 
disappears for the zigzag junction when the G?2~wave sym- 
metry is chosen. We therefore speculate that the origin 
of the 'bulge' is the response of the on-site order param- 
eter to the lack (or in the case of d\ + jc?2-wave, partial 
lack) of order parameter overlap and/or the suppression 
in supercurrent that effectively follows from this. 



E. LDOS and the existence of ZES 

Finally we discuss the LDOS for different junctions. 
Figure [D shows the LDOS for a clean zigzag SNS junc- 
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tion with J = (a) and J ^ with a di-wave symmetry 
in N (b) and for a clean armchair SNS junction with 
J = (c) and J ^ with a d\ + ic?2-wave symmetry in 
N (d). First it is clear that the abrupt drop in the ef- 
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FIG. 9: LDOS for zigzag junction with L = 10 unit cells and 
J — (a), J = 4 eV (b) and for armchair junction with 
L = 104 unit cells and J = (c), J = 4 eV (d), all with 
/x(N) = 0.7 eV. Zigzag junction (b) has di-wave symmetry in 
N which has nodes in the gap function, armchair junction (d) 
has di -Mob-wave symmetry which is fully gapped. Black color 
corresponds to 2 states/eV/unitcell. Vertical lines indicates 
the interfaces. 

fective chemical potential parameter at the SN interface 
gives an abrupt drop even in the electron filling of the 
bands. This is seen in the lighter color in the N region 
which corresponds to a lower DOS and is the result of 
the N region being closer to the zero DOS Dirac point. 
Second we see that for a relatively short junction (a) the 
energy gap is not depleted inside N when J — whereas 
for a longer junction (c) it is to a large extend fully de- 
pleted due to the length of the junction. When J ^ 
we see in (b) the appearance of the nodal quasiparticle 
spectrum characteristic of the di-wave superconducting 
state whereas in (d) we see a large energy develop due to 
the fully gapped d\ + id,2-wave state. 

There is some tendency for interfacial states in the gap 
in Fig. [D^d) but the DOS is still very low and no local- 
ized states are actually formed. This however raises the 
question of the existence of zero energy states (ZESs) at 



surfaces or interfaces. For high-T c superconductors it was 
shown relatively early on that midgap surface states will 
exist at certain surfaces^ Later on this was extended to 
also include the formation of ZESs at various interfaces 
between superconductors with unconventional order pa- 
rameters (see Ref. [2{| and references therein). The sur- 
face state can be viewed as a bound state formed in a N 
layer on top of the superconductor in the limit where the 
thickness of N goes to zero. Bound states are formed in 
this layer due to the retro-reflectivity of the Andreev re- 
flection which causes an electron traveling in the N layer 
to form a closed trajectory. For conventional supercon- 
ductors these states coincides with the de Gennes-Saint- 
James bound states/22, However, for unconventional or- 
der parameters the two Andreev reflections involved in 
this process will take place via different order parame- 
ters. In order for a surface state to form at zero energy 
it is only necessary that A(0) and A(ir — 9), where 9 is 
the angle of incidence for an electron quasiparticle at the 
surface, have different signs^ In our case the necessary 
sign change in the order parameter is true for the e?2-wave 
symmetry at a zigzag surface and the c?i-wave symmetry 
at an armchair surface where the latter case appears nat- 
urally for a clean interface, as seen in Fig. [5] Despite this 
we were not able to detect any trace of a ZES at the arm- 
chair surface. We attribute the absence of a ZES to the 
local symmetry character changes in the surface. Within 
the first few unit cells of the armchair surface the charac- 
ter of the order parameter has both strong s- and G?2~wave 
components, which neither give rise to a ZES. In fact, 
the only time we were able to see the signature of a ZES 
was for a zigzag surface when we forced the symmetry 
to be G?2-wave in the whole slab. A similar enforcement 
of di-wave symmetry in a whole armchair slab leads to a 
depleted pair potential at the surface, and consequently 
also a lack of a ZES. It might be worth emphasizing that 
we are here dealing with rather heavily doped graphene 
samples so we do not expect the rather peculiar specu- 
lar Andreev reflection possible near the Dirac points to 
possibly play a role here. 

ZESs in a SS' junction have a similar origin to the 
surface states. In fact, in the tunneling limit the exis- 
tence of a ZES is determined independently at each of 
the two superconductor surfaces. In the opposite, fully 
transparent limit, a sign change in the order parameters 
has to be present across the junction instead. In the re- 
gion of intermediate transparency a linear combination of 
these two conditions applies. 29 For a square lattice with 
d 2 .2_j / 2-wave superconductor junctions treated within the 
TB-BdG formalism, ZESs were reported for several in- 
terfaces including the {110} interface: 3 - However, it was 
also found that when using a lattice model the appear- 
ance of a ZES is sensitive to Friedel oscillations in the 
wave function. These can cause destructive interference 
between different surface lattice sites, leading to the dis- 
appearance of the ZES. For our graphene junctions with 
conventional s-wave on one side and d-wave on the other 
side, the same criterion as for a surface for seeing a ZES 
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applies, independent of the transparency of the interface. 
For the junctions at hand, only the rough interface e^- 
wave symmetry zigzag junction would qualify but, as in 
the surface case, the character of the order parameter is 
very mixed at the interface, and we do not see any sign 
of a ZES. 

Finally, let us comment on the case of the time-reversal 
symmetry breaking d\ + id2-wave state. This state will 
not have any ZES, but it can still have subgap localized 
states, which can be viewed as a splitting of the zero 
energy peak^ However, apart from the small DOS at the 
interface in Fig. OJd) we have not been able to see any 
signs of subgap states. The Andreev reflection process 
described above to give rise to the ZES at the surface of 
the superconductor will in fact in this case produce a net 
current along the interface. This is, however, beyond the 
scope of this paper. 



IV. CONCLUSIONS 

We have investigated the possibility of enhancing the 
intrinsic nearest neighbor spin-singlet (SB) correlations, 
present in all p7r-bonded planar organic molecules, in 
graphene by constructing a graphene SNS Josephson 
junction with conventional s-wave superconducting con- 
tacts. For strong enough SB coupling constant J (and 
doping level) a d\ + ie^-wave superconducting state de- 
velops in the bulk. However, we have found that in a 
finite-sized junction, as in a finite-sized slab, a larger J 
than in the bulk is needed to achieve superconductiv- 
ity. In addition, the SB superconducting state will still 
have d-wave symmetry despite the nonzero coupling be- 
tween the SB extended s-wave state to the extrinsic s- 
wave state induced from the superconducting contacts. 
Using a realistic estimation of J, it would be necessary 
to have a very high doping level in the N region in order 
to develop a d-wave state. While such high doping levels 
cannot as of present be achieved by traditional gating of 
the graphene, there might be a possibility of using chem- 
ical doping, such as is assumed for the S regions in this 
work. However, care has to be taken to not change any 
other physical characteristics of the graphene other than 
the Fermi level. An interesting side remark is that, for a 
dangling bond short armchair slab, it might be possible 
to see a SB superconducting state for lower doping levels 
than in the bulk thanks to a strong singlet formation on 
the surface bond. 

We have also found that interface effects are very im- 
portant for the choice of d-wave symmetry. For even the 
longest clean zigzag junctions we could model a di-wave 
symmetry is favored whereas for the clean armchair junc- 
tion <i2-wave symmetry is favored in shorter junctions 
but the di + id2-wa,ve state is reached for sufficiently long 
junctions. Also for the clean junctions, the on-site s-wave 
pair amplitude in S and the d-wave in N do not influence 
each other and the situation can to good approximation 
be described as a 'bulk-meets-bulk' junction. For the d\- 



wave in a zigzag junction or the c^-wave in an armchair 
junction there is also no coupling or overlap between the 
d-wave in N and the s-wave in S, inhibiting any Joseph- 
son coupling in the junction. In fact, the supercurrent 
decreases when the d-wave develops in N as quasiparticle 
tunneling and Andreev reflection are suppressed due to 
the nodal spectrum developing in the quasiparticle DOS 
in N for the d\- and d2-wave. The d± + zd2-wave has a 
finite overlap of one of its components with the S region s- 
wave state but the quasiparticle DOS in N is fully gapped 
and, when both effects are combined, we have found that 
the net effect is still a decrease in the current. We have 
also reported on the development of a 'bulge' in the phase 
of the on-site pair amplitude when a d-wave develops in 
N. This structure is pronounced but will not influence 
the current vs phase relationship in the junction and it is 
only present when (part of) the d-wave state lacks cou- 
pling with the on-site s-wave. For the clean zigzag or 
armchair interfaces this 'bulge' feature is, except for a 
< 10% decrease in current, the only sign of a d-wave in 
N in the usual Josephson junction characteristics. The 
structure is however pronounced and might be possible 
to detect with a phase sensitive measurement inside N. 

In order to facilitate coupling between the two pair 
amplitudes we have also investigated interfaces with 
smoother doping profiles and atomic scale interface 
roughness and been able to see an increased on-site pair 
amplitude, and thus higher supercurrent, when the d- 
wave state develops in N. However, in order to drasti- 
cally increase the current, the opposite d-wave symmetry 
than the one found naturally in clean junctions has to 
be favored. We have shown that this is possible when 
including interface roughness and the current has been 
shown to increase with 50% in one such case. This would 
be a clear signal for a developing superconducting state 
due to SB correlations. 

Because of computational constraints we have only in- 
vestigated the two simplest interfaces, zigzag and arm- 
chair, but in a general experimental junction it is likely 
that a chiral and also nonsmooth interface is present. 
Due to the near degeneracy in energy of the two different 
d-wave solutions and their complex combinations, it is 
impossible to say, without detailed knowledge about the 
interface, which symmetry will be favored in a general 
junction and as seen here, the choice of symmetry will 
have a big influence on the properties of the junction. 
However, as an experimental junction will at least expe- 
rience some doping leakage into N and not be perfectly 
clean on the atomic level, it will, with high enough dop- 
ing in the N region to induce a SB superconducting state, 
show an increased Josephson current due to an increased 
on-site pairing potential, and possibly even a significant 
increase due to a particular choice of d-wave symmetry 
in N. 

Despite the unconventional order parameters and the 
theoretical existence of zero energy states (ZESs) at a 
surface or interface of such superconductors, most junc- 
tions and slabs we have investigated do not show any 
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signs of ZESs in the gap. We attribute this to the del- 
icate symmetry character mixing at surfaces and inter- 
faces. Therefore, it is unlikely that one can detect the 
presence of a d-w&ve state in an experimental system by 
searching for a zero energy peak in the LDOS or the zero 
bias conductance peak it is known to create. 

Finally, let us point out the possible benefit of a SNS 
junction with <i-wave contacts, made by depositing, for 
example, a high-T c material on graphene. As has been 
clear from this study, the biggest problem with conven- 
tional contacts is the zero coupling between the s-wave 
and the d-wave states in the bulk. For rough interfaces, 
coupling can to some extent be obtained but the funda- 
mental symmetry difference is still a big obstacle. While 
it remains to be studied, we believe that ci-wave contacts 
could potentially greatly enhance the effect of the in- 



trinsic SB correlations. A corresponding study has been 
made on the square lattice with d-wave contacts and 
a N region with d-wave correlations but T c < T and 
pronounced effects in proximity effect and current were 
seen^ For graphene there would also be the additional 
issue of having multiple d-wavc symmetries which com- 
plicate the picture, but if any effects could be seen, these 
would experimentally prove the existence of SB correla- 
tions in graphene. 
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